Radial distribution function of penetrable sphere fluids to second order in density 
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The simplest bounded potential is that of penetrable spheres, which takes a positive finite value 
e if the two spheres are overlapped, being otherwise. In this paper we derive the cavity function 
to second order in density and the fourth virial coefficient as functions of T* = fesT/e (where fcs 
is the Boltzmann constant and T is the temperature) for penetrable sphere fluids. The expressions 
are exact, except for the function represented by an elementary diagram inside the core, which is 
approximated by a polynomial form in excellent agreement with accurate results obtained by Monte 
Carlo integration. Comparison with the hypernetted-chain (IfNC) and Percus-Yevick (PY) theories 
shows that the latter is better than the former for T* ^ 1 only. However, even at zero temperature 
(hard sphere limit), the PY solution is not accurate inside the overlapping region, where no practical 
cancelation of the neglected diagrams takes place. The exact fourth virial coefficient is positive for 
r* < 0.73, reaches a minimum negative value at T* « 1.1, and then goes to zero from below as 
1/r** for high temperatures. These features are captured qualitatively, but not quantitatively, by 
the HNC and PY predictions. In addition, in both theories the compressibility route is the best one 
for T* < 0.7, while the virial route is preferable if T* > 0.7. 

PACS numbers: 61.20.Gy, 61.20.Ne, 05.20.Jj, 05.70.Ce 



I. INTRODUCTION 

Ultrasoft and bounded potentials represent useful 
models to characterize the effective two-body interaction 
in some colloidal systems, such as star or chain polymers 
in good solvents P, i, i, i, i, H, 0, H- The simplest 
bounded potential is that of so-called penetrable spheres 
(PS), which is defined as 



e, r < cr, 
0, r > (T, 



(1.1) 



where e > 0. This interaction potential was suggested 
by Marquest and Witten Q as a simple theoretical 
approach to the explanation of the experimentally ob- 
served crystallization of copolymer mesophases and it 
has been since then the subject of a number of stud- 
ies 0, E [H E, El, Q El, El, E3, El, El M, EH- 

The classical integral equation theories, in particular the 
Percus-Yevick (PY) and the hypernetted-chain (HNC) 
approximations, do not describe satisfactorily well the 
structure of the PS fluid, especially inside the overlap- 
ping region for low temperatures. Thus, the PS model 
can be used as a stringent benchmark to test alterna- 
tive theories jH, El, El ES [SI. From that point of 
view, the derivation of exact properties provides an in- 
valuable tool. The exact structural and thermodynamic 
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properties of the PS fluid in the high-temperature limit 
T* = UbT / e —f oo (where ks is the Boltzmann constant 
and T is the temperature) are known for any density pa^ , 
including the high-density regime pa^ ~ T* [ll| . On the 
other hand, the corresponding properties in the comple- 
mentary low-density limit for any temperature has not 
been addressed, to the best of our knowledge, except in 
the one-dimensional case i21]. 

The aim of this paper is to derive the exact expressions 
for the radial distribution function g(r) and, equivalently, 
the cavity function y{r) of PS fluids to second order in 
density. To that end we will exploit the fact that the PS 
Mayer function is proportional to the hard sphere (HS) 
Mayer function. This implies that the diagrams to be 
evaluated are the same as in the case of HS, except that 
now each diagram is affected by a temperature-dependent 
factor. 

In the next Section we present some definitions and 
basic equations. The density expansion of y{r) to second 
order is worked out in Sec. IIIIl where the HS functions 
derived by Nijboer and van Hove [11] outside the core 
r > a are complemented by their extensions in the over- 
lapping region (r < cr). However, we have not been able 
to derive the rigorously exact expression for r < cr of 
the function x(r) represented by the only elementary di- 
agram. Instead, the exact values of x(0), x'(0), x(f), 
x'(ct), x"(f), x"'(f), and Jq dr r"^ xi''') s^-re obtained in 
Sec. IIVI With these constraints, we have constructed a 
polynomial approximation of x(r) for r < a which yields 
results indistinguishable from those obtained by Monte 
Carlo (MC) integration with six significant figures. The 
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exact fourth virial coefRcient is also derived in Sec. IIVI 
The exact results are compared with the HNC and PY 
predictions in Sec. |Vl It is seen that the latter is gen- 
erally preferable at low temperatures, while the former 
is more accurate at high temperatures. The paper ends 
with the conclusion section. 



The series expansions of the cavity function and the 
compressibility factor in powers of density read 

oo / fi \ " 

y(rh,T*) = l + 5]y„(r|T*) - ^\ (2.9) 

n=l ^""^ 



II. DEFINITIONS AND BASIC EQUATIONS 

We consider in this paper a fluid of particles interacting 
via the pairwise potential . Henceforth we take a = 
1 as the length unit. Let us introduce the cavity (or 
background) function 



2/(rh,r*) = e*W/'^--^5(H'7,T*), 



(2.1) 



where g{r\rj^T*) is the radial distribution function, 77 = 
(7r/6)/9 being the packing fraction. Equation (|2.ip implies 
that 

5(r|ry, T*) = y{r\T^, T*) ~ xyir\ri, r*)e(l - r), (2.2) 

where Q{r) is the Heaviside step function and we have 
called 
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(2.3) 



The parameter x represents the probability of rejecting 
an overlap of two particles in a MC move. The thermody- 
namic quantities can be expressed in terms of g(r\ri,T*) 
or y{r\r], T*) Particularized to the PS model, 

the compressibility factor Z = p/ pksT is given by the 
virial equation of state as 



Z(77,r*) ^l + Aijxv{l\ij,T*). 



(2.4) 



The (dimensionless) isothermal compressibility K 
ksT {dp /dp) J, is 



1 + 2477 



drr^ [y{r\r],T*) ~ 1] 



X I drr y{r\r],T*) 



(2.5) 



Finally, the internal energy per particle can be written 



uir,,T*) 



-T* + Ut]{1 - x) [ drr^y{r\T],T*) 
2 Jo 



(2.6) 

These three quantities are thermodynamically connected 
by the relations 



djyZ) 

drj 



(2.7) 



Z{v,T*) = l 



(2.10) 



In Eq. dmni), 6„(r*) is the (reduced) nth virial coef- 
ficient. The quantities {6„(T*)} can be obtained from 
the functions {yn{r\T*)} through the virial route, Eq. 
(12. 4p . the compressibility route, Eq. (|2.5p . or the en- 
ergy route, Eq. (|2.6p . In order to distinguish the re- 
sults derived through each route, we will use the no- 
tation bl{T*), b^„{T*), 6^(r*), respectively. Of course, 
bl{T*) = b^{T*) = bf,{T*) if the exact cavity function is 
employed. 

Insertion of the expansion (|2.9|) into Eqs. (|2.4p and 
(1^ yields (for n > 2) 



(2.11) 



bl{T*) = l2{n-l)^-^ dx, drrV-2(r|rr). 

(2.12) 

In Eq. (|2.12p use has been made of Eq. (|2.8p and of the 
ideal gas condition limT*-too bn{T*) = 0. In the case of 
the compressibility route, insertion of Eq. (|2.9p into Eq. 
(|2.5p and use of the relation (|2.7p leads to the recursive 
formula 



n-l 



where Ki{T*) 
for n> 2. 



(2.13) 



m— 1 

-Sx and 
24(^ 

\ 7T 



/ drr^j/„_i(r|r* 



-X drr y„-i(r|r*) 
Jq 



(2.14) 



III. CAVITY FUNCTION TO SECOND ORDER 
IN DENSITY 

The virial coefficients yn{r\T*) are represented by di- 
agrams m, [2^. In particular. 



d{u/e) 
' dr] 



(1 



) — 

dx 



(2.8) 



yiir\T*) 




(3.1) 
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y2{r\T*) 



(3.2) 

Here, the open circles represent root points separated 
by a distance r, the fiUed circles represent field points 
to be integrated out, and each bond represents a Mayer 
function 

/(r|T*) = e-^W/*=«'^ - 1. (3.3) 
Thus, for instance, 



rfr3/(ri3|T*)/(r23|T*), (3.4) 





= / ^1-3 y dr4/(ri3|T*)/(r34|T*) 

x/(r24|T*)/(ri4|T*), (3.5) 



where — \vi — Yj\ and ri2 = r. 

Equations (I3.1|) - (I3.5|) hold for any interaction poten- 
tial. In the special case of PS, the Mayer function be- 
comes 



where 



/(r|r*)=x/Hs(0, 



/Hs(0 = -e(i-r) 



(3.6) 



(3.7) 



is the Mayer function of HS. Therefore, the spatial de- 
pendence of each one of the diagrams contributing to the 
virial expansion (|2.9p is exactly the same as for HS. The 
only difference is that each diagram is now multiplied 
by the temperature-dependent parameter x raised to a 
power equal to the number of bonds in that particular 
diagram. As a consequence, Eqs. p.ip and (|3.2p become 



yi{r\T*) = x^7(r), 



(3.8) 
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y2{r\T*) = xV(r) + 2x^(0 + y7'(0 + yX(0- (3-9) 

Here, "f(r) is represented by the diagram on the right- 
hand side of Eq. (|3.ip . except that now each bond corre- 
sponds to a Mayer function /hs ■ Analogously, the func- 
tions ip(r), ^iid x(^) are represented by the first, 



second, and fourth diagram, respectively, on the right- 
hand side of Eq. (|3.2p , with /hs for each bond. The ex- 
pressions of these functions for r > 1 are known [2^ . [2^ . 
The region r > 1 is the physically relevant one in the case 
of HS. However, the overlapping region r < 1 is essential 
in the case of PS, since g{r\ri,T*) ^ for r < 1, except in 
the zero-temperature limit (where the PS model reduces 
to the HS one). Therefore, it is necessary to extend the 
knowledge of 7(r), (/?(r), ijj{r), and xi''') to the domain 
< r < 1. 

Given a radial function F{r) we define its Fourier 
transform as 



dre^'-'Fir) 



An 

T 



dr r sm(kr)F{r) 



It is easy to realize that j{k) = [/Hs(fc)]^, where 

~ , , , k cos k — sin k 

/Hs(fc) = 47r- 



k^ 

Inverse Fourier transform simply yields 

7(r) = ^(2-r)2(r + 4)e(2-r). 



(3.10) 



(3.11) 



(3.12) 



This implies that the function "f{r) for < r < 1 is just 
the analytical continuation of its expression for 1 < r < 
2. Next, note that (p{k) — [fasik]]'^, so that 



(^(r) = (^A(r)e(l -r) + (^B(r)e(3 - r) 



(3.13) 



with 



^sir) = -^^{r - 3)\r' + 12r2 + 27r - 6). (3.15) 

Therefore, (p{r) — (pA{r) + 'fB{i') for Q < r < \, while 
Lp{r) = ^b{t) for I < r < 3. In the case of V'('')j one 
has ^{k) = /Hs(fc)7*(fc), where 7*(r) = 7(0/hs(0- As 
a consequence. 



V'(r) = ?AA(r)e(l - r) + VjB(r)e(2 - r) (3.16) 



with 



V'a(0 = ~-ipA{r), 



(3.17) 



1 



^jbM = ^(r-2)2(r-'^+4r''-51r^-10r2+479r-81). 

36 35r 

(3.18) 



Equations (|3.15p and (|3.18p coincide with those derived 
in Ref. ^22] by a different method. On the other hand, 
the functions (pAij") and ipAij"), which are needed to get 
Lp{r < 1) and ip{r < 1), respectively, were not considered 
in Refs. [13] and [1^. Near the origin. 



7(r) 



6 



(8-6r)+0(r2), 



(3.19) 
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^{r) = -—m + 0{r^), 



ip{r) 



36 



(30- 15r) + 0(r2). 



(3.20) 



(3.21) 



Equations (|3.13[) - (|3.18[) show that ip{r) has a fourth- 
order discontinuity at r = 1 and at r = 3, while ip{r) 
has a fourth-order discontinuity at r = 1 and a second- 
order discontinuity at r = 2. 

Now we turn to the much more involved function xC*"); 
represented by the elementary diagram at the end of the 
right-hand side of Eq. (|3.2p . Let us decompose it in a 
form similar to Eqs. (|3.13|) and (|3.16p . 

Xir) ^ XAir)e{l~r)+XBir)e{VS-r)-^^{r). (3.22) 

The exact expression for Xsif) was obtained by Nijboer 
and van Hove 22]. It reads 



XB{r) 



/3r2 
I 280 



41 
420 



V3(4-r2) 



3-r2 


-( 


23 

15" - 


36 
35^ 


3r6 






2r 


560 " 


15 


+ y 


^15 



cos 



9 
35^ 

15 



r-3 



3r6 
560 



r 

15 



9 
35^ 



V3(4-r2) 



(3.23) 



We have not been able to obtain an analytic expression 
for x(^) ill the interval < r < 1. By working with 
bipolar coordinates, it is possible to express the derivative 
x'{r < 1) as a sum of 13 triple integrals, but only two 
of them seem to be analytically solvable. Therefore, we 
have resorted to numerical evaluation of < 1) by 
the MC method [1^ and to a very accurate polynomial 
approximation. In order to construct the latter, some 
constraints on the exact xi''' < 1) Siie derived in the next 
Section. 



IV. CONSTRAINTS ON x{r)- POLYNOMIAL 
APPROXIMATION 



From Eqs. ((3?22)) and ((3?23)) one can get 
' 36 I 2 4 / ' 



TT^ 1 / 3A7bf^ 7219 256V2 \ 
36 51 3 6 ~ J 



(4.1) 



(4.2) 



X"(l) ^ - ^ + '-^] , (4.3) 



36 153 \ 3 6 



TT / 



^'"(1) ^ J_ |^ 946bP _ 16597 _ 4832^2 ^^ ^ ^^^^^ 



36 153 \ 3 3 

In the above equations, 

2707 438%/2-4131sec-i3 



6f = 



70 



707r 



18.3648 (4.5) 



is the exact value of the fourth virial coefficient for HS. 
Next, note that 

X(0) - J dvj{r)fus{r)^~Y I dr r^r - 2f {r + 4) 



- -36^0- 



(4.6) 



The same result is obtained from the following zero- 
separation theorem for HS [27| : 

lnyHs(0|??) = 477?;Hs(lh) + 4 /" rfryiyHs(l|?7i), (4-7) 

where UHsiflv) — 1™t*^o y(''|?7, 2^*) is the cavity func- 
tion for HS. As a further constraint on the unknown 
function x(r) for r < 1, let us consider the alternative 
zero-separation theorem 



yHs(oi^) 
yHs(0|?7) 



-6tohs(1|'7)- 



(4.8) 



This implies that limT-^o ^2(017"*) = -(7rV36)63. From 
Eqs. dSll) and ((37T9| - ([37ST|l one then has 



In this Section we derive some constraints on xi''') for 
r < 1. First, we take into account that x(r) and its first 
three derivatives must be continuous at r = 1. We are 
not aware of a formal proof of this statement, but it is 
strongly supported by the following two arguments: (i) 
both ip{r) and ipir) have a fourth-order discontinuity at 
r — 1, even though a diagonal bond is added when going 
from the diagram representing ip{r) to that representing 
'0(r); (ii) in the one-dimensional case, the three functions 
(p{r), ip{r), and x(^) have the same type of s ing ularity at 
r = 1, namely a second-order discontinuity [2l[. 



X'(0) = ^30. 



(4.9) 



As a consequence of Eqs. ((3?T9l) - ((3?2T1) . (|46| . and (|49|) . 
the form of y2 near the origin is 

y2ir\T*) = — a;M-30 + 2a;(46- 39r)- 15x^(1 -r)] 
36 



(4.10) 



Let us apply now the condition of thermodynamic con- 
sistency for the fourth virial coefficient 64 (T*). Taking 



5 



into account that yoir\T*) = 1, yi{l\T*) = x'^jil) 
(57r/12)a;2, and 



36" 



544 6347 f bf^ 57 



Eq. (P1T|) yields 



64 (T*) 



(1 



4352 - 1995a; 



70 



(4.11) 



(4.12) 



(4.13) 



The same results for &2 and 63 are obtained through the 
energy route, Eq. (|TT^ . As for 64, Eq. (|TT^ yields 



b4T*) = ^9x^ [ dr 
^ Jq 



ip{r) + -xip{r) 



+ -x-f'^{r) + -x'^xir) 
5 3 



(4.14) 



Since the functions ^{r), '4'{r), and 7^(r) arc known for 
r < 1, the integrals involving them can be performed. 
Thus, equating the right-hand sides of Eqs. (|4.13|) and 
(I4.14P one gets 



6HS 57 



drr^Y(r) = — I — 

^ 36 V 3 6 



In turn, this condition implies that 



IxH). 



(4.15) 



12752 
256 -—X - 



6347 



x^ + {2bf^ - 57) x^ 



35 35 

(4.16) 

where the coefficients Kn arc defined by Eq. (|2.14p . Use 
of Ki = -8x, K2 = 2x^(32 - 152:), and (jil^ in Eq. 
(1^1^ leads again to Eqs. (jil^ and 28]. 

Although the exact analytic expression of Xa{'>'), and 
hence of xi'"') for r < 1 is not known, we have derived 
in this Section a number of constraints. The value of 
x(r) and its first three derivatives at r = 1 are given by 
Eqs. (|iT|) - ((i:i)) . On the other hand, Eqs. (jiTB)) and (g^ 
give x{t) and x' {f) at the origin. Finally, the integral of 
r^x(r) in the interval < J' < 1 is determined by Eq. 
(j4.15p . Since there are seven constraints we can approx- 
imate xi'"') for r < 1 by a polynomial of sixth degree: 



36 



"0 



+ ai{r- 1) -f a2(r- 1)^ 



+a3(r - 1)3 + (r ~ 1)^ (/3o + /3ir + (i^r^)] . 

(4.17) 



In this equation, the constants = (36/7r^)x(l) 
ai = {m/^^)x'{l), a2 = {m/^^)x"{l)/2, and as = 



(36/7r^)x"'(l)/6 are obtained from Eqs. 
From Eqs. (|4.6p and (|4.9p one gets 



(l4T])-(l4a 



1 / 43096f s 
^° - 459 



129317 10784V2\ 

6 TT / ' 



^ 1 / 554&fs 8663 1120\/2\ 
Pi = I 



27 \ 3 3 7 

Finally, application of (|4.15p yields 

1 / 38036^^3 134713 920\/2 



6 



12 



(4.19) 



(4.20) 



The second derivative at the origin is 
„ _ TT^ 10 / 550696JS 981592 22232\/2^ 



36 



2.07929. 



(4.21) 



In contrast, the exact result is (see the Appendix) 



- — 2.07608. (4.22) 

Ok) 



Therefore, Xpoiy(0)/x"(0) ^ 1.00155. This gives an idea 
of the extreme accuracy of Xpoiy('')- In fact, we have 
evaluated numerically x(^) by MC integration with 6 sig- 
nificant figures and have found that Xpoiy('") agrees with 
x{r) within the error bars (see Table |T|. One could ex- 
ploit the exact knowledge of x"(0), Eq. (|4.22p . to propose 
an approximation presumably even more accurate than 
Eq. (|4.17p , but this does not seem to be necessary in view 
of Table H 



V. COMPARISON WITH THE HNC AND PY 
THEORIES 

Once we have obtained the exact temperature- 
dependence of the function y2{'r) and the associated 
fourth virial coefficient 64 for PS, it is worthwhile compar- 
ing these two quantities with the predictions provided by 
the two classical integral equation theories, namely the 
HNC and PY theories. 



A. Cavity function to second order, y2{r) 

In the HNC theory, the elementary diagrams are ne- 
glected at any order in density [25]. To second order 
in density, the only elementary diagram is the last one 
given in Eq. (|3.2p . Therefore, the function y2{r) is ap- 
proximated by 



(5.1) 



6 



TABLE I: Values of -(36/7r^)x(r) in the interval < r < 1 
as obtained numerically by MC integration and as given by 
the polynomial approximation (|4.17p . The number enclosed 
between parentheses in the second column indicates the 95%- 
confidence error. 



T 


MC 


Eq. (|4.17l) 




0.00 


29.9994(6) 


30 




0.05 


28.5029(5) 


28.5031 




0.10 


27.0146(5) 


27.0144 




0.15 


25.5367(4) 


25.5369 




0.20 


24.0735(4) 


24.0736 




0.25 


22.6277(4) 


22.6275 




0.30 


21.2015(3) 


21.2017 




0.35 


19.7987(3) 


19.7991 




0.40 


18.4228(3) 


18.4228 




0.45 


17.0756(3) 


17.0758 




0.50 


15.7614(2) 


15.7612 




0.55 


14.4822(2) 


14.4820 




0.60 


13.2414(2) 


13.2413 




0.65 


12.0420(2) 


12.0421 




0.70 


10.88745(13) 


10.88740 




0.75 


9.78032(12) 


9.78037 




0.80 


8.72413(14) 


8.72404 




0.85 


7.72152(10) 


7.72151 




0.90 


6.77582(8) 


6.77586 




0.95 


5.89018(7) 


5.89019 




1.00 


5.06759(5) 


5.06762 





In the PY approximation, apart from the elementary di- 
agrams, a subset of the remaining diagrams is also ne- 
glected. In particular, the PY expression for y2{r) only 
retains the two first diagrams in Eq. (|3.2|) . so that 



(5.2) 



Figure [T] compares the exact function y2 (r) with the 
HNC and PY approximations at T* = (hard spheres), 
T* = 1, and T* = 2. Both theories agree very well 
with the exact 2/2 (^) for r > 1.5 but discrepancies are 
apparent for shorter distances, especially inside the core 
(r < 1). Although restricted to low densities. Fig. [1] 
clearly illustrates some of the general features found at 
finite densities [H, [13] : the HNC overestimates the pen- 
etrability effect, while the PY approximation underesti- 
mates it. The former property is a consequence of the 
neglect of (x^ /2)x{r), which is a negative definite quan- 
tity. This is only partially compensated by the PY ne- 
glect of (a:^/2)7^(r), since 7^(r) > |x('')| for r < 1 and, 
moreover, > x^. While the PY theory tends to be 
better at lower temperatures (i.e., when the overlapping 
of particles is hindered and the system is close to that 
of HS), the HNC is preferable at higher temperatures. 
If we characterize the quality of each approximation by 
the separation of the corresponding contact value ^2(1) 



from the exact result, it turns out that the temperature 
beyond which the HNC approximation becomes better 
than the PY approximation is T* ~ 1.04. This is similar 
to the behavior found in the onc-dimcnsional case 12 ll. 



15 
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r 

FIG. 1: (Color online) Plot of the function y2{r) at T* = 
(top panel), T* = 1 (middle panel), and T* = 2 (bottom 
panel). The solid lines are the exact results, the dashed lines 
are the HNC predictions, and the dotted lines are the PY 
predictions. 



B. Fourth virial coefficient 

The knowledge of yf^'^{r\T*) and y2^{r\T*) allows 
one to get the associated expressions for the fourth virial 



7 



TABLE II: Fourth virial coefficient 64 (T*) and other related quantities as given exactly and by the HNC and PY approximations 
through the virial (v), energy (e), and compressibility (c) routes. 



Theory b^T') MO) TfT &4|„,,„ 

Exact x*[bf^x^ - {1 - x){4352 - 1995x)/70] b^^ 0.7250 1.1027 -1.4803 

HNC,v/e a;'*(6347a;- 4352)/70 ^ 0.8641 1.2574 -1.1258 

HNC,c a;'' (31735a; - 261 12)/420 ^ 0.5778 0.9314 -2.3345 

PY,v 16x'*(171a; - 136)/35 16 0.6304 0.9888 -2.0378 

PY,e 2a;''(6347a; - 5440)/175 ^ 0.5140 0.8641 -2.7485 

PY,c x''(6347s - 4352)/105 19 0.8641 1.2574 -0.7505 



coefficient b4{T*). As discussed in Section [TTl there 
are three alternative routes [cf. Eqs. (|2.1ip - (|2.13p ] and 
there is no reason to expect internal consistency among 
them, unless the exact y2ir) is used. The expressions 
for bliT*), bl{T*), and bi(T*) in the HNC and PY 
approximations are given in Table |TT1 where, for com- 
pleteness, the exact expression, Eq. (|4.13p . is also in- 
cluded. It is known ^29*| that the HNC integral equation 
provides thermodynamically consistent results through 
the virial and energy routes, regardless of the potential 
considered. This explains the fact that &J'^^'"(T*) — 



On the other hand, the PY integral equa- 



tion yields three different predictions, i.e., b^^'^{T*) ^ 
(T*). It is interesting to note that 
(T*) = ^xdbl''-^'iT*)/dx. 



lPY,, 



'4 (T*) ^ 

HNC,v/ef 
4 

In the limit T* one recovers the known results 

for HS, namely bT'''"^"{0) = f = 28.5, ^^^''(O) 



2 "4 



^ ~ 13.3881, 64 '"(0) = 16, and 64 '"(0) = 19. Al- 
though the energy route is ill defined for strict HS, taking 
the zero-temperature limit on the PS model yields well 
defined values [s^. In that way, one finds 6r'"(0) = 
~ 10.3657, which is a rather poor value reflect- 
ing the inaccuracy at any temperature of 2/2^(7") for 
r < 1. In the opposite high-temperature limit, one 
has limT-^oo64(T*)/a;4 = -2176/35. This exact value 
is retained by all the approximations, except by the 
compressibility route in the PY theory, which yields 
limT-^00 bl^'''{T*)/x'^ = -4352/105, i.e., 2/3 of the ex- 
act result. 

While 62 (T*) and 63 (T*) are positive definite quanti- 
ties, this is not the case of 64 (T*). The latter quantity 
changes sign at a certain "Boyle- like" temperature Tq. 
In addition, bi[T*) presents a (negative) minimum value 
^4 1 mill at a temperature Tmin > ^0 • The numerical values 
of Tq*, T*jin, and bi\^^^^ are displayed in Table [TTl From 
Eq. (|2.12p one can see that the temperature Tj^j,, asso- 
ciated with 5|(T*) is the temperature across which the 
integral dr r'^y2{r\T*)^ and hence the third-order term 
in the density expansion of the internal energy, changes 
from positive to negative. 

The temperature dependence of the fourth virial coef- 
ficient is shown in Fig. [21 where, apart from the exact 
curves, the two HNC approximations and the three PY 



-q' 10 



-20 




FIG. 2: (Color online) Plot of the fourth virial coefficient 
b4,{T*) (top panel) and of the scaled fourth virial coefficient 



-l/T* 

exactly and by the HNC and PY approximations. 



hi{T*)/x*' (bottom panel), where a; = 1 — e 



as given 



approximations are included. The best approximation 

PY c 

up to T* ~ 0.71 is provided by 64 ' . In the intermedi- 
ate range 0.71 ^ T* < 1.04, however, 64^''' presents the 
best agreement. Finally, for T* ^ 1.04 the best perfor- 
mance corresponds to 6^'^^'"^'^. Within the PY theory, 
the energy route is never better than the virial route 
but becomes preferable to the compressibility route for 
T* > 1.22. In the case of the HNC theory, the compress- 
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ibility route is better than the virial/energy routes for 
T* < 0.73 only. 



VI. CONCLUSION 

In this paper we have considered a three-dimensional 
fluid of particles interacting via the PS interaction (jl.ip . 
This potential encompasses the ideal gas in the high- 
temperature limit (T* — > oo or a; ^ 0) and the HS fluid 
in the low-temperature limit (T* ^ or a; — *■ 1). How- 
ever, at finite temperature the problem becomes much 
more difficult. Even the one-dimensional case is not ex- 
actly solvable 21] since there is no a priori limitation to 
the number of particles that can interact simultaneously 
with a given particle. 

The diagrams which appear in the density expansions 
for the PS fluids are exactly the same as those appearing 
for HS fluids, except that each diagram needs to be mul- 
tiplied by the temperature-dependent parameter x raised 
to the number of bonds. By exploiting this fact, we have 
obtained the cavity function through second order in den- 
sity and the equation of state through the fourth virial 
coefficient. In order to obtain y2{r\T*), we have needed 
to extend to r < 1 the functions evaluated by Nijboer 
and van Hove [22] for r > 1. Nevertheless, the possible 
analytical evaluation of the elementary-diagram function 
x(r) for r < 1 seems to be a formidable task. Thus, 
we have resorted in that case to two complementary ap- 
proaches: (i) a numerical computation by MC integration 
with an error bar of the order of 0.001% and (ii) a sixth- 
degree polynomial approximation constructed by enforc- 
ing seven exact constraints. Both methods show such 
an excellent mutual agreement that the results obtained 
from the polynomial approximation can be considered as 
exact from a practical point of view. 

The results obtained here for y2{r\T*) and 64 (T*) have 
been compared with those corresponding to the two clas- 
sical integral equation theories, namely the HNC and PY 
theories. It is known that the PY theory is much better 
than the HNC one for HS fluids, so that one could have 
expected a similar situation for PS fluids, at least at low 
temperatures. Our results show that this is indeed the 
case, provided that T* < 1. However, even at very low 
temperatures (including the HS limit T* 0), the PY 
solution strongly underestimates the cavity function in 
the overlapping region. This reflects the fact that the 
fortunate practical cancelation (in the case of HS) of the 
diagrams neglected by the PY equation does not apply 
for r < 1. In this respect, it is interesting to note that 
the widely extended belief that the PY theory becomes 
exact in the special case of one-dimensional hard rods is 
only correct for r > 1 [2l| . 

When comparing the exact fourth virial coefficient 
with the HNC and PY theories one has to take into 
account their thermodynamic inconsistency, yielding 
two HNC predictions (virial/energy and compressibil- 
ity routes) and three PY predictions (virial, energy, and 



compressibility routes) . All these predictions capture the 
non-monotonic behavior of b^{T*). In both theories, the 
compressibility route is the best one for T* < 0.7, while 
the virial route is preferable if T* > 0.7. As in the case 
of the structural functions, the equation of state is better 
described by the HNC equation than by the PY equation 
for high enough temperatures (T* ^ 1). 

It is obvious that access to non-trivial exact informa- 
tion on the structural and thermodynamic properties of 
fluids, even if restricted to special cases, is of paramount 
importance. From that point of view, we hope that the 
results reported in this paper can contribute to an ad- 
vancement on our knowledge of the behavior of systems 
of particles interacting through bounded potentials. 



Acknowledgments 

One of the authors (Al.M.) is grateful to the Junta de 
Extrcmadura for supporting his stay at the University 
of Extremadura in the period October-December 2005, 
when this work was started. His research has been par- 
tially supported by the Ministry of Education, Youth, 
and Sports of the Czech Republic under the project No. 
LC 512 and by the Grant Agency of the Czech Repub- 
lic under project No. 203/06/P432. The research of 
the other author (A.S.) has been supported by the Min- 
isterio de Educacion y Ciencia (Spain) through Grant 
No. FIS2004-01399 (partially financed by FEDER funds) 
and by the European Community's Human Potential 
Programme under contract No. HPRN-CT-2002-00307, 
DYGLAGEMEM. 



APPENDIX: EVALUATION OF x"(0) 

The function x(r) is represented by the elementary di- 
agram displayed at the end of the right-hand side of Eq. 
(jX^ . Thus, 

X{r2)^ j drs jdrif{r:i)f{ri)f{r2z)f{r2A)f{r3A), 

(A.l) 

where here f{r) — fusif) = —0(1 — r). Now we dif- 
ferentiate with respect to r2 and take into account the 
mathematical property 

=(5(r23-l)-^ — = (5(r23-l) ■ A. 2) 

dr2 dr2 r2 

The result is 

X{r2) = 2jdr3j dr4/(r3)/(r4)/(r24)/(r34) 

X(5(r23-l)cos023, (A.3) 

where ^23 is the polar angle of the vector r23 and the z 
axis is assumed to point in the direction of r2. Making 
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the change of variables r23, r4 r24, Eq. (|A.3p 

becomes 



X'(r2) = 2j drsj dr4/(r23)/(r24)/(r4)/(r34) 

x5{r3 - 1)00363. (A. 4) 



Note that r^g 



1 — 2r2cos03, so that a necessary 



(but not sufficient) condition for /(r'23) ^ is 613 < 7r/2. 
Therefore, in the Hmit r2 ^ 0, one has 



;^'(0) = 47r7(l) / d6'3sin6'3Cos6l3 
Jo 



36 



30, (A.5) 



in agreement with Eq. (|4.9p . 

Now we differentiate again with respect to r2 to get 



X"(^2)-x'l'(r-2)+X2(^2), 



(A.6) 



where 



X'i'(r2) ^ 2j dr^J dr4/(r24)/(r4)/(r34) 

x5(r3 - l)cos6'3(5(r23 - 1) cos 6*23, (A.7) 



In Eq. (|A.8|) , since = r| + 1 — 2r2 cos ^24, a necessary 
condition for /(r4) ^ is 6*24 < 7r/2. Now, setting r2 = 
and taking into account that cos6'24 — cos 6*4, Eq. 
(|A.8|) becomes 



xi'(O) = -2 / d03 / 
Jo Jo 

X I c?(cos 6*4 ) cos 6*4 / (r34 ) 



4 / (i(cos6'3)cos03 

JO 



(A.13) 

where now r^^ = 2[1 — cos 63 cos ^4 — sin 63, sin ^^4 cos((/)3 — 
^4)]. The changes z = COS63, z' = — cos^?4, and = 
^3 — 04 lead to 



X2(0) = -Stt [ dz [ dz'zz' 

Jo Jo Jo 



x9 cos( 



1 + 2zz' 



2^{l-z^){l-z'^) 



(A.14) 



X'2'(r2) = 2j dv3j dVif{r^3)f{rA)f{r3i) 

xS{r3 - 1) cos6'3(5(r24 - 1) cos6'24- (A. 8) 

Let us first consider Xi ('")• Note that cos 6*23 = r2— cos 6*3, 
where it has been taken into account that = 1. Now, 
using the property 

S{h{x)) = \h'{xo)\-^S{x-xo), (A.9) 

where h(x) is a function that vanishes at x = a;o, we have 

S{r23 - 1) = r2 i,5(cos03 - ^2/2). (A.IO) 

Thus, 



Xi(r2) = 2r2 



3 / dr4/(r24)/(r4)/(r34), (A.ll) 



- 1 — r4[r2 cos 6*4 + ■y/4 — r2 sin 64 cos( 



where r^^ = 

04)], 03 and 04 being azimuthal angles. At the origin 
one simply has 



It can be easily seen that 1 + 2zz' < 2y (1 — z^){l — z' ) 

if and only if 2^ + z'^ + zz' < 3/4. This requires that 
< z < v/372 and < z' < (^3(1 - z^) - z)/2. Conse- 
quently, 



nV3/2 A^Z{\~Z-^)^7.)I2 

X'2'(0) = -Stt / dz \ dz'zz' 
Jo Jo 

1 + 2zz' 



X cos 



2ja-z^)ii-z'') 



(A.15) 



The result of the integral is 



(A.16) 



X'i'(0) = 0. 



(A.12) 
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